#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string.h>
#include "header.h"

/*	
	Multiple traits analysis ver.05

	This program can treat ordinary multiple-trait linear additive genetic models and structural equation models (SEMs)

	For SEMs, the structure among traits and statistical models can fixed specified with "Structure file".
	For example, 

		0 0 0
		2 0 0
		0 1 0

	Columns indicate the upstream (affecting) traits and rows indicate the downstream (affected) traits.
	The numbers indicate the models.
	In this example, trait 1 affects trait 2 with model 2, and trait 2 affects trait 3 with model 1.

	This program recognizes the relationship between trait 1 and 2 as the "first", trait 2 and 3 as the "second" relationships.
	For example, output files for MCMC samples of lambda, are named as "Lambda1" and "Lambda2" respectively.

	Current models:

		1: lambda is shared among all records (homogeneity).
		The prior distribution is 

			lamdba ~ Normal (0, sigmar)

		sigmar is set to 10000 arbitrary.

		2: a Dirichlet process prior is used as a prior distribution of lambda of each record.
		The Dirichlet process is

			lambda ~ Dirichlet process ( F(lambda), alpha )

		where F(Lambda) is the base distribution and alpha is the concentration parameter.
		F() is the normal distribution with parameters specified by the user.

		3: B-spline 

			Lambda[y] = B.0(y)*P[0] + B.1(y)*P[1] + ... + B.Nsp-1(y)*P[Nsp-1]

		where B.0 ~ B.Ncp-1 are the base functions defined by deBoorCox, Ps are the coefficients, and Nsp is the number of splines.
		When SplinePrior==1, the prior distributions of the spline coefficients are

			P[0] ~ Normal (0, BecomeVague*sigmaPi)
			P[i] ~ Normal (P[i-1], sigmaPi) i>=1

		When SplinePrior==2, the prior distributions of the coefficients are

			P[0] ~ Normal (0, BecomeVague*sigmaPi)
			P[1] ~ Normal (0, BecomeVague*sigmaPi)
			P[i] ~ Normal (2*P[i-1]-P[i-2], sigmaPi) i>=2

		where BecomeVague = 1000 (default). The prior of sigmaPi is scaled inverse chi square.

	The linear model can fixed written as

		Y = LY + XB + ZU + E

	where Y is observations, L is the recursive or simultaneous effects, B is fixed effects or covariates, X and Z are the design matrices, 
	U is the additive effects and E is residuals.

	Z is common for all traits.
	X is unique to each trait.

	The relationship matrix for additive effects should fixed given as an inverted matrix.

	Input files are as follows.

		Observations file	: observations should fixed sorted by the order of trait and record. K (number of traits) * N (number of records) records included.
		Va file				: includes K * K elements of the scaling matrix for the prior inverted Wishart distribution of genetic (co)variances.
		Ve file				: includes K * K elements of the scaling matrix for the prior inverted Wishart distribution of environmental (co)variances.
								If "0", residual covariance is not assumed or residual correlation is fixed to the values given by ResidualCor file
		ResidualCor file	: defines correlation of residuals among traits. Non-informative prior is assigned.
								If "0", residuals are assumed to be independent among traits. Non-informative prior is assigned.
		parameter file		: defines parameters.
		Structure file		: includes K * K elements. This defines the sturcture among traits.
		X file				: design (covariates) matrix X. X includes K matrices for trait 1 to K, which are placed horizontally in the file.
								When the total number of fixed effects are F, X is N * F matrix.
		Z file				: design matrix Z. Traits share the single matrix. When the number of additive effects is P, Z is N * P matrix.
								If "0", an identity matrix is used. In this case, N is equal to P.
		Relationship file	: includes P * P elements of the inverted relationship matrix.

*/

/* Random seed. Time is used when -1 */
#define  seed -1

/* Objects about this algorithm */ 
#define  Algoname  "mta05"					/* name of this algorithm */

int main (int argc, char *argv[])
{
	/* Values defined by the parameter file */
	char	AnalysisName[101]="";				/* Analysis name. Output files are named as "Algoname_AnalysisName_..." */
	char	Observationsfile[101]="";			/* Observations file name */
	char	Relationshipfile[101]="";			/* Inverted relationship matrix file name. Ignored if Vafile = 9*/
	char	Vafile[101]="";						/* Va file name */
	char	Vefile[101]="";						/* Ve file name */
	char	ResidualCorfile[101]="";			/* File defines residual correlation. */
	char	Xfile[101]="";						/* X file name. Ignored when sum(Ft)=0 */
	char	Zfile[101]="";						/* Z file name. */
	int		Leaveoneout;						/* If 0, models are fitted. If n (>0), leave one out validation starts with the nth individual. Now under construction*/
	double	Missingvalue;						/* Missing records in the observation file */
	int		K;									/* Number of traits. When increase K, check the value of KKlimit in the header file. */
	int		P;									/* Number of additive effects */
	int		N;									/* Number of records per trait */
	double	va;									/* Degree of freedom for the prior inverted Wishart distribution of genetic (co)variances */
	double	ve;									/* Degree of freedom for the prior inverted Wishart distribution of environmentl (co)variances. Ignored if Vefile==0 */
	double	*Nu;								/* Nu of the prior inverse chi square distribution of residual variance. Used when CovR=0 or 1. */
	double	*S2;								/* S2 of the prior inverse chi square distribution of residual variance. Used when CovR=0 or 1. */
	int     Burnin;								/* Length of burnin */
	int     Ite;								/* Length of iterations after burnin */
	int     Thi;								/* Frequency of sampling */
	int		*Ft;								/* Number of fixed or covariates effects for each trait */
	int		SplinePrior;						/* Define the prior distributions of the spline coefficients */
	int		AddIntercept;						/* If 1, the first basis spline is replaced by the intercept. Used in model 3 */
	char	Structurefile[101]="";				/* Structure file name. If "0", reduced (ordinary multi-trait models) are applied */

	/* In addition, the parameter file should contain the additional information (DiricletPrior and NspandOrd) at the last row(s) when model 2 or 3 are selected .

	Model2
		DirichletPrior Meanbd Sigmabd Alpha
	where Meanbd and Sigmabd denote the mean and the variance of the base distribution, respectively, and Alpha is the concentration parameter of DP.

	Model3
		NspandOrd Nsp Ord Nupi S2pi
	where Nsp is the number of splines, Ord is the order of the base functions, and Nupi and S2pi are the parameters of scaled inverse chi square 
	distributions of sigmaPi.

	*/

	/* Seed used */
	int		usedseed;

	/* Inputs and outputs */
	char    Input[101]="", Output[101]="", Labeltemp[51]="";
	double	*Invertedmatrix, *Outputmatrix;
	FILE    *fp, *fp2;

	/* For repeat statement */
	int		ii, jj, mm;
	int		trait, record, additive, fixed;

	/* Generic use */
	double	temp, temp2;

	/* Observations */
	Ystruct	*Y;

	/* Assumption on covariance. 0: residuals are independent among traits, 1: residual correlation is fixed, 2: residual covariance is inferred */
	int		CovR;

	/* Lambda */
		/* Impute recursive or simultaneous effects, or not (SEMs or ordinary multi-trait models)*/
		int		SEM;

		/* network structure, number of recursive or simultanesous effects */
		int		*Structure, M;

		/* Structure of Lambda */
		Lambdastruct *L;

		/* mean and variance of the prior distribution of Lambda, which is used when model 1.*/
		double	Lambda0 = 0.0, Lambdasigmafixed = 10000.0;

	/* Location parameters */
		/* Structure for fixed effects */
		Lstruct	*X, *Z;

		/* Total number of fixed or covariates effects, cumulative number of effects */
		int		F, *cumFt;

		/* Whether the additive effects are assumed or not*/
		int		Polygenic;

		/* use an identity matrix as a design matrix (Z) */
		int		Identitymatrix;

	/* Variance components */
		/*inverted relationship matrix */
		double	*iA;

		/* Va matrix, Ve matrix (prior of covariance structure)*/
		double	*Va, *Ve;

		/* Residual correlation */
		double	*ReCor;

	/* For sampling */
		/* Sampled parameters */
		double	*SamplediG, *SamplediR, *SampledLikelihood, *MeanR, *iMeanR, MeanLoglike;

	/* Products of varables */
	int		KN, KK, PP12, KK12, PK;

	/* Information criterion */
	double	LoglikeFromPostMean, PDIC, DIC, WAIC, PWAIC2, LPPD;
	double	*Dummy, *Products, *Eri;

	/* For sampling parameters */
	int		Ns, Nswb, Nsab;

	/* Statistics of observations */
	double	*Mean, *Sd;

	/* Open the parameter file */
	if ( ( fp = fopen( argv[1], "r" ) ) == NULL )
	{
	   	printf("can not open the parameter file (%s)\n", argv[1]);
		return (0);
	}

	/* Read the parameter file */
	fscanf (fp, "%s%s", Labeltemp, AnalysisName);
	fscanf (fp, "%s%s", Labeltemp, Observationsfile);
	fscanf (fp, "%s%s", Labeltemp, Relationshipfile);
	fscanf (fp, "%s%s", Labeltemp, Vafile);

		if (strcmp(Vafile, "9")==0){ Polygenic=0; } else { Polygenic=1; }

	fscanf (fp, "%s%s", Labeltemp, Vefile);
	fscanf (fp, "%s%s", Labeltemp, ResidualCorfile);

		if (strcmp(Vefile, "0")==0){
			if (strcmp(ResidualCorfile, "0")==0){
				CovR = 0; /* residuals are independent among traits */
			}else{
				CovR = 1; /* residual correlation is defined by ResiudalCorfile */
			}
		}else{
			CovR = 2; /* residual covariance is inferred. The prior is defined by Vefile */
		}

	fscanf (fp, "%s%s", Labeltemp, Xfile);
	fscanf (fp, "%s%s", Labeltemp, Zfile);

		if (strcmp(Zfile,"0")==0) {Identitymatrix=1;} else {Identitymatrix=0;}

	fscanf (fp, "%s%d", Labeltemp, &Leaveoneout);
	fscanf (fp, "%s%lf",Labeltemp, &Missingvalue);
	fscanf (fp, "%s%d", Labeltemp, &K);

		KK = K*K;
		if (K == 1) CovR = 0; /* when K==1, CovR should be fixed to 0 */

	fscanf (fp, "%s%d", Labeltemp, &P);
	fscanf (fp, "%s%d", Labeltemp, &N);

		if (strcmp(Zfile,"0")==0)
			if ( P != N ){printf ("P should fixed equal to N, when Zfile is 0 (identity matrix is used)\n"); return(0);}

	fscanf (fp, "%s%lf", Labeltemp, &va);
	fscanf (fp, "%s%lf", Labeltemp, &ve);

		Nu = (double*) calloc(K, sizeof(double));
		S2 = (double*) calloc(K, sizeof(double));

	fscanf (fp, "%s", Labeltemp);
	for(trait=0; trait<K; trait++) fscanf (fp, "%lf",Nu+trait);
	fscanf (fp, "%s", Labeltemp);
	for(trait=0; trait<K; trait++) fscanf (fp, "%lf",S2+trait);

	fscanf (fp, "%s%d", Labeltemp, &Burnin);
	fscanf (fp, "%s%d", Labeltemp, &Ite);
	fscanf (fp, "%s%d", Labeltemp, &Thi);
	fscanf (fp, "%s", Labeltemp);

		Ft = (int*) calloc ( sizeof(int), K );		if(Ft==NULL) alert();
		cumFt = (int*) calloc ( sizeof(int), K );	if(cumFt==NULL) alert();
		F = 0;
		for (trait=0; trait<K; trait++)
		{
			Ft [trait] = -10;
			fscanf (fp, "%d", Ft+trait);	
			F += Ft[trait];
			if (trait>0) cumFt[trait] = cumFt[trait-1] + Ft[trait-1];
		}

		if (Ft[K-1]==-10) { printf ("the length of Ft is lower than K (%d)\n", K); return(0);}

		ii=-10;	fscanf (fp, "%d", &ii);
		if (ii!=-10) { printf ("the length of Ft is higer than K (%d)\n", K); return(0);}

	fscanf (fp, "%s%d", Labeltemp, &SplinePrior);
	fscanf (fp, "%s%d", Labeltemp, &AddIntercept);
	fscanf (fp, "%s%s", Labeltemp, Structurefile);

		if ( strcmp(Structurefile,"0")==0 ) { SEM=0; } else { SEM=1;}
		if ( SEM==1 && K==1) { printf ("SEM is chosen although K == 1\n"); return(0);}

		/* Open the Structure file and read additional arguments if model 2 or 3 is specified */
		if (SEM)
		{
			sprintf ( Input, "%s", Structurefile );
			if ( ( fp2 = fopen( Input, "r" ) ) == NULL )
			{
				printf("can not open the structure file\n" );
				return (0);
			}
			Structure = (int*) calloc ( KK, sizeof ( int ) ); if (Structure == NULL) alert();

			M = 0;
			for (ii=0; ii<K; ii++ )
			{
				for (jj=0; jj<K; jj++ )
				{
					fscanf ( fp2, "%d", Structure+ii*K+jj );
					if ( Structure [ii*K + jj] < 0 || Structure [ii*K + jj] > 3) { printf ( "Values in the structure file should fixed 0, 1 or 2\n"); return(0);}

					if (Structure[ii*K+jj] > 0)
						M ++;
				}
			}
			fclose (fp2);

			L = (Lambdastruct*) calloc (M, sizeof(Lambdastruct));		if(L==NULL) alert();
			mm=0;
			for (ii=0; ii<K; ii++)
			{
				for (jj=0; jj<K; jj++)
				{
					if (Structure[ii*K+jj] > 0)
					{
						L[mm].update = Structure [ii*K+jj];
						L[mm].down = ii;
						L[mm].up = jj;
						mm ++;
					}
				}
			}

			/* Meanbd, Sigmabd, Alpha */
			fscanf (fp, "%s", Labeltemp);
			for (mm=0; mm<M; mm++)
				switch(L[mm].update){
				case 2:
					fscanf (fp, "%lf %lf %lf", &(L[mm].meanlambda),&(L[mm].lambdasigma),&(L[mm].alpha));
					if (L[mm].lambdasigma<=0.0){printf ( "Sigmabd should fixed >0.0 (now %f)\n", L[mm].lambdasigma ); return(0);}
					if (L[mm].alpha<=0.0){printf ( "Alpha should fixed >0.0 (now %f)\n", L[mm].alpha ); return(0);}
					break;
				case 3:
					fscanf (fp, "%d %d %lf %lf", &(L[mm].Nsp),&(L[mm].Ord), &(L[mm].Nupi), &(L[mm].S2pi));
					if ((L[mm].Nsp-L[mm].Ord)<1){printf ( "Nsp (%d) - Ord (%d) should satisfy >=1\n", L[mm].Nsp, L[mm].Ord ); return(0);}
					break;
			}
	}

	fclose (fp);

	/* Products of parameters */
	KN = K*N;
	PK = P*K;
	PP12 = P*(P+1)/2;
	KK12 = K*(K+1)/2;

    /* Number of samples */
    Ns = (Burnin + Ite) / Thi;
	Nswb = Burnin / Thi;
	Nsab = Ns - Nswb;

	/* Open the observation file */
	sprintf ( Input, "%s", Observationsfile );

	if ( ( fp = fopen( Input, "r" ) ) == NULL )
	{
	   	printf("can not open the observations file\n" );
		return (0);
	}

	Y = (Ystruct*) calloc ( K, sizeof ( Ystruct ) ); if (Y == NULL) alert();
	for (trait=0; trait<K; trait++)
	{
		Y[trait].observations = (double*) calloc (N, sizeof(double)); if(Y[trait].observations==NULL) alert();
		for ( record=0; record<N; record++ )
			if ( fscanf ( fp, "%lf", Y[trait].observations + record ) == -1) 
			{
				printf ( "The number of elements in the observations file is fewer than KN (%d)\n", KN );
				putchar ( '\n' ); return (0);
			}
	}
	if (fscanf ( fp, "%lf", &temp ) > 0) 
	{
    	printf ( "The number of elements in the observations file is larger than KN (%d)\n", KN);
    	putchar ( '\n' ); return (0);
	}
	fclose (fp);

	/* Open the relationship matrix file */
	if ( Polygenic )
	{
		sprintf ( Input, "%s", Relationshipfile );

		if ( ( fp = fopen( Input, "r" ) ) == NULL )
		{
	   		printf("can not open the relationship matrix file\n" );
			return (0);
		}
		iA = (double*) calloc ( PP12, sizeof ( double ) ); if (iA == NULL) alert();
		mm = 0;
		for (ii=0; ii<P; ii++ )
			for (jj=0; jj<P; jj++ )
			{
				if( fscanf ( fp, "%lf", &temp ) == -1)
				{
					printf ( "The number of elements in the relation matrix file is fewer than P*P (%d)\n", P*P );
					putchar ( '\n' ); return (0);
				}

				if (jj>=ii){ iA[mm] = temp; mm++;}
			}

		if ( fscanf ( fp, "%lf", &temp ) > 0 )
		{
    		printf ( "The number of elements in the relation matrix file is larger than P*P (%d)\n", P*P);
    		putchar ( '\n' ); return (0);
		}
		fclose (fp);
	}

	/* Open the Vafile */
	if ( Polygenic )
	{
		Va =(double*) calloc ( KK, sizeof (double) );

		if (K > 1)
		{
			if (Vafile[0] == '0' )
			{
				for (ii=0; ii<K; ii++)
					Va [ ii*K + ii ] = 1.0;
			}
			else
			{
				sprintf ( Input, "%s", Vafile );
				if ( ( fp = fopen( Input, "r" ) ) == NULL )
				{
					printf("can not open the Vafile\n" );
					return (0);
				}

				for (ii=0; ii<KK; ii++ )
				{
					if (fscanf ( fp, "%lf", Va+ii ) == -1)
					{
						printf ( "The number of elements in the Vafile is fewer than K*K (%d)\n", KK );
						putchar ( '\n' ); return (0);
					}
				}

				if ( fscanf ( fp, "%lf", &temp ) > 0 )
				{
    				printf ( "The number of elements in the Vafile is larger than K*K (%d)\n", KK);
    				putchar ( '\n' ); return (0);
				}
				fclose (fp);
			}
		} /*if (K > 1)*/
	}

	/* Open the Vefile */
	if(CovR == 2){
		Ve = (double*) calloc ( KK, sizeof ( double) );		if(Ve == NULL) alert();
		if (Vefile[0] == '0' )
		{
			for (ii=0; ii<K; ii++)
				Ve [ ii*K + ii ] = 1.0;			
		}
		else
		{
			sprintf ( Input, "%s", Vefile );
			if ( ( fp = fopen( Input, "r" ) ) == NULL )
			{
	   			printf("can not open the Vefile\n" );
				return (0);
			}
			for (ii=0; ii<KK; ii++ )
			{
				if (fscanf ( fp, "%lf", Ve+ii ) == -1)
				{
					printf ( "The number of elements in the Vefile is fewer than K*K (%d)\n", KK );
					putchar ( '\n' ); return (0);
				}
			}

			if ( fscanf ( fp, "%lf", &temp ) > 0 )
			{
    			printf ( "The number of elements in the Vefile is larger than K*K (%d)\n", KK);
    			putchar ( '\n' ); return (0);
			}
			fclose (fp);
		}
	}else{ /* When CovR == 0 or CovR == 1 */

		if(CovR == 1){
			/* Open ResidualCorfile */
			ReCor = (double*) calloc ( KK, sizeof ( double) );		if(ReCor == NULL) alert();
			sprintf ( Input, "%s", ResidualCorfile );
			if ( ( fp = fopen( Input, "r" ) ) == NULL )
			{
	   			printf("can not open the ResidualCorfile\n" );
				return (0);
			}
			for (ii=0; ii<KK; ii++ )
			{
				if (fscanf ( fp, "%lf", ReCor+ii ) == -1)
				{
					printf ( "The number of elements in the ResidualCorfile is fewer than K*K (%d)\n", KK );
					putchar ( '\n' ); return (0);
				}
			}

			if ( fscanf ( fp, "%lf", &temp ) > 0 )
			{
    			printf ( "The number of elements in the ResidualCorfile is larger than K*K (%d)\n", KK);
    			putchar ( '\n' ); return (0);
			}
			fclose (fp);
		}
	}

	/* Open X file */
	if (F>0)
	{
		sprintf ( Input, "%s", Xfile );
		if ( ( fp = fopen( Input, "r" ) ) == NULL )
		{
	   		printf("can not open the Xfile\n" );
			return (0);
		}

		/* Read X file and convert it to a proper form */
		X = (Lstruct*) calloc (F, sizeof(Lstruct));														if(X==NULL) alert();
		for (fixed=0; fixed<F; fixed++) { X[fixed].covariates = (double*) calloc (N, sizeof(double));	if(X[fixed].covariates==NULL) alert(); }
		for (ii=0; ii<N; ii++)
		{
			for (fixed=0; fixed<F; fixed++)
			{
				if(fscanf (fp, "%lf", X[fixed].covariates+ii)==-1)
				{
					printf ( "The number of elements in the X file is fewer than F*N (%d)\n", F*N );
					putchar ( '\n' ); return (0);
				}
			}
		}

		if (fscanf ( fp, "%lf", &temp ) > 0 )
		{
    		printf ( "The number of elements in the X file is larger than F*N (%d)\n", F*N);
    		putchar ( '\n' ); return (0);
		}
		fclose (fp);
	}

	/* Open Z file */
	if (Polygenic)
	{
		Z = (Lstruct*) calloc ( PK, sizeof ( Lstruct ) );															if (Z == NULL) alert();
		for (additive=0; additive<P; additive++){ Z[additive].covariates = (double*) calloc (N, sizeof(double));	if (Z[additive].covariates==NULL) alert();}
		/* covariates are loaded to Z for the first trait, because covariates are common for all traits */

		if (Zfile[0] == '0') 
		{
			for (additive=0; additive<P; additive++ )
				Z[additive].covariates[additive] = 1.0;
		}
		else
		{
			sprintf ( Input, "%s", Zfile );
			if ( ( fp = fopen( Input, "r" ) ) == NULL )
			{
	   			printf("can not open the Zfile\n" );
				return (0);
			}

			/* read Z file */
			for (record=0; record<N; record++)
				for (additive=0; additive<P; additive++)
					if (fscanf ( fp, "%lf", Z[additive].covariates + record) == -1)
					{
						printf ( "The number of elements in the Z file is fewer than P*N (%d)\n", P*N );
						putchar ( '\n' ); return (0);
					}

			fscanf ( fp, "%lf", &temp );
			if ( fscanf ( fp, "%lf", &temp ) > 0 )
			{
    			printf ( "The number of elements in the Z file is larger than P*N (%d)\n", P*N);
    			putchar ( '\n' ); return (0);
			}
			fclose (fp);			
		}
	}

/*----modify from here for LOO ------------------*/
	/* Initialize the seed */
	usedseed = Randomize ( seed );
	putchar ( '\n' );
	printf ( "Random seed: %d\n", usedseed );
	putchar ( '\n' );
	ix = ( rand()%30000 ) + 1;
	iy = ( rand()%30000 ) + 1;
	iz = ( rand()%30000 ) + 1;

	/* Output the parameters */
	sprintf ( Output, "%s_%s_parameters.txt", Algoname, AnalysisName );

	fp = fopen ( Output, "w" );
	fprintf ( fp, "Observationsfile : %s\n", Observationsfile );
	fprintf ( fp, "Relationshipfile : ");

	if ( Polygenic ) { fprintf ( fp, "%s\n", Relationshipfile);} else { fprintf( fp, "Ignored\n");}

	fprintf ( fp, "Vafile           : " );

	if ( Polygenic==0 ) { fprintf ( fp, "9 (no polygenic effects)\n"); } else { if (K==1) { fprintf (fp, "Ignored (noninformative inverse-chi square)\n");} else { if ( Vafile[0] == '0' ) { fprintf ( fp, "0 (identity matrix)\n"); } else { fprintf ( fp, "%s\n", Vafile);} } }

	fprintf ( fp, "Vefile           : %s\n", Vefile );
	fprintf ( fp, "ResidualCorfile  : %s\n", ResidualCorfile );
	fprintf ( fp, "CovR             : %d", CovR);
	
	switch(CovR){
	case 0:
		fprintf ( fp, " (independent residuals)\n");
		break;
	case 1:
		fprintf ( fp, " (correlation is fixed)\n");
		break;
	case 2:
		fprintf ( fp, " (covariance is inferred)\n");
		break;
	}

	fprintf ( fp, "Xfile            : " );

	if ( F>0) { fprintf ( fp, "%s\n", Xfile); } else { fprintf ( fp, "Ignored\n"); }

	fprintf ( fp, "Zfile            : ");

	if (Polygenic==0) { fprintf ( fp, "Ignored\n"); } else { if ( Zfile[0] == '0') { fprintf ( fp, "0 (identity matrix)\n"); } else { fprintf ( fp, "%s\n", Zfile); }}

	fprintf ( fp, "Missingvalue     : %f\n", Missingvalue );
	fprintf ( fp, "K                : %d\n", K );
	fprintf ( fp, "P                : %d\n", P );
	fprintf ( fp, "N                : %d\n", N );
	fprintf ( fp, "va               : ");
	if ( K==1 ) { fprintf (fp, "Ignored\n"); } else { fprintf (fp, "%f\n", va);}

	fprintf ( fp, "ve               : ");
	if ( CovR==2 ) { fprintf (fp, "%f\n", ve); } else { fprintf (fp, "Ignored\n");}

	fprintf ( fp, "Nu               : ");
	if (CovR==2) fprintf (fp, "Ignored\n");
	for(trait=0; trait<K; trait++) fprintf (fp, "%f ", Nu[trait]);
	fprintf (fp, "\n");

	fprintf ( fp, "S2               : ");
	if (CovR==2) fprintf (fp, "Ignored\n");
	for(trait=0; trait<K; trait++) fprintf (fp, "%f ", S2[trait]);
	fprintf (fp, "\n");

	fprintf ( fp, "Burnin           : %d\n", Burnin );
	fprintf ( fp, "Ite              : %d\n", Ite );
	fprintf ( fp, "Thi              : %d\n", Thi );
	fprintf ( fp, "Ft               : " );
	for (trait=0; trait<K; trait++)
		fprintf ( fp, "%d ", Ft[trait]);
	fprintf ( fp, "\n");
	fprintf ( fp, "SplinePrior      : %d\n", SplinePrior );
	fprintf ( fp, "AddIntercept     : %d\n", AddIntercept );
	if (SEM)
	{
		fprintf ( fp, "Structurefile    : %s\n", Structurefile );
		for(mm=0; mm<M; mm++)
		{
			fprintf (fp, "Lambda%d Model%d", mm+1,L[mm].update);
			switch(L[mm].update)
			{
			case 1:
				fprintf (fp, "\n");
				break;
			case 2:
				fprintf (fp, " %f %f %f\n", L[mm].meanlambda, L[mm].lambdasigma, L[mm].alpha);
				break;
			case 3:
				fprintf (fp, " %d %d %f %f\n", L[mm].Nsp, L[mm].Ord, L[mm].Nupi, L[mm].S2pi);
				break;
			}
		}
	}
	else
	{
		fprintf ( fp, "An ordinary multiple-trait model is applied\n");
	}

	fprintf ( fp, "Seed             : %d\n", usedseed );
	fclose (fp);

	/* For sampling */
	if (Polygenic)
	{
		for (trait=0; trait<K; trait++)
			for (additive=0; additive<P; additive++)
			{
				Z[trait*P + additive].sampledeffect = (double*) calloc (Ns, sizeof (double));		if(Z[additive].sampledeffect == NULL) alert();
			}
		SamplediG = (double*) calloc ( KK12 * Ns, sizeof (double) );								if(SamplediG==NULL) alert();
	}

	if (CovR==0) 
	{ 
		SamplediR = (double*) calloc ( K * Ns, sizeof (double));				if(SamplediR==NULL) alert();
	}
	else
	{
		SamplediR = (double*) calloc ( KK12 * Ns, sizeof (double) );			if(SamplediR==NULL) alert();
	}

	if (F>0) 
	{
		for (fixed=0; fixed<F; fixed++)
		{
			X[fixed].sampledeffect = (double*) calloc (Ns, sizeof(double));		if(X[fixed].sampledeffect==NULL) alert();
		}
	}

	for (trait=0; trait<K; trait++)
	{
		Y[trait].meanerrors = (double*) calloc (N, sizeof(double));				if(Y[trait].meanerrors==NULL) alert();
	}

	if (SEM)
	{
		for (mm=0; mm<M; mm++)
		{
			switch (L[mm].update) {			
			case 1:
				L[mm].sampledlambda = (double*) calloc ( Ns, sizeof(double) );				if(L[mm].sampledlambda==NULL) alert();
				break;
			case 2:
				L[mm].sampledlambda = (double*) calloc ( N*Ns, sizeof(double));				if(L[mm].sampledlambda==NULL) alert();
				L[mm].sampledclass = (int*) calloc ( N*Ns, sizeof(int));					if(L[mm].sampledclass==NULL) alert();
				break;
			case 3:
				L[mm].sampledPi = (double*) calloc ( Ns*L[mm].Nsp, sizeof(double) );		if(L[mm].sampledPi==NULL) alert();
				break;
			}
		}
	}

	/* likelihood calculation */
	SampledLikelihood = (double*) calloc ( Ns*N, sizeof(double));					if(SampledLikelihood==NULL) alert();

	/* mean and sd of observations */
	Mean = (double*) calloc (K, sizeof(double));
	Sd = (double*) calloc (K, sizeof(double));

	/* Start MCMC sampling */
	mta05 (Missingvalue, K, P, N, va, ve, Burnin, Ite, Thi, CovR, Ft, F, cumFt, Y, L, X, Z, SEM, M, Lambda0, Lambdasigmafixed,
	Polygenic, iA, Va, Ve, ReCor, Nu, S2, Identitymatrix, Mean, Sd, SamplediG, SamplediR, SampledLikelihood, SplinePrior, AddIntercept);

	/* Output */
	/* SamplediG */
	Invertedmatrix = (double*) calloc ( KK, sizeof(double) );		if (Invertedmatrix==NULL) alert();
	Outputmatrix = (double*) calloc ( KK, sizeof(double) );			if (Outputmatrix==NULL) alert();

	if (Polygenic)
	{
		sprintf ( Output, "%s_%s_G.txt", Algoname, AnalysisName );
		fp = fopen ( Output, "w" );
		for ( ii=0; ii<Ns; ii++ )
		{
			if (K==1)
			{
				fprintf (fp, "%f\n", 1.0/SamplediG[ii]);// * Sd[0] * Sd[0]);
			}
			else
			{
				for ( trait=0; trait<K; trait++ )
				{
					for ( jj=0; jj<K; jj++ )
					{
						if (jj>=trait)
						{
							Invertedmatrix [trait*K + jj] = SamplediG [ ii * KK12 + trait*K - (trait-1)*trait/2 - trait + jj];
						}
						else
						{
							Invertedmatrix [trait*K + jj] = SamplediG [ ii * KK12 + jj*K - (jj-1)*jj/2 - jj + trait];
						}
					}
				}
				GaussJordan ( Invertedmatrix, Outputmatrix, K);
				for (trait=0; trait<K; trait++)
					for (jj=0; jj<K; jj++)
						fprintf (fp, "%f ", Outputmatrix[trait*K+jj]);// * Sd[trait] * Sd[jj]);
				fprintf ( fp, "\n" );
			}
		}
		fclose (fp);
	}

	/* SamplediR */
	MeanR = (double*) calloc ( KK, sizeof (double) );
	iMeanR = (double*) calloc ( KK, sizeof (double) );
	sprintf ( Output, "%s_%s_R.txt", Algoname, AnalysisName );
	fp = fopen ( Output, "w" );
	for ( ii=0; ii<Ns; ii++ )
	{
		if (CovR==0 )
		{
			for ( jj=0; jj<K; jj++)
			{
				temp = 1.0/SamplediR [ii*K + jj];
				fprintf ( fp, "%f ", temp);// * Sd[jj] * Sd[jj]);
				if (ii > Nswb )
					MeanR [jj * K + jj] += temp;
			}
		}
		else
		{
			for ( trait=0; trait<K; trait++ )
			{
				for ( jj=0; jj<K; jj++ )
				{
					if (jj>=trait)
					{
						Invertedmatrix [trait*K + jj] = SamplediR [ ii * KK12 + trait*K - (trait-1)*trait/2 - trait + jj];
					}
					else
					{
						Invertedmatrix [trait*K + jj] = SamplediR [ ii * KK12 + jj*K - (jj-1)*jj/2 - jj + trait];
					}
				}
			}
			GaussJordan ( Invertedmatrix, Outputmatrix, K);
			for (trait=0; trait<K; trait++)
				for (jj=0; jj<K; jj++)
				{
					fprintf (fp, "%f ", Outputmatrix[trait*K+jj]);// * Sd[trait] * Sd[jj]);
					if (ii > Nswb )
						MeanR [trait*K+jj] += Outputmatrix[trait*K+jj];
				}
		}
		fprintf ( fp, "\n" );
	}
	fclose (fp);

	for (jj=0; jj<KK; jj++)
		MeanR [jj] /= (double)Nsab;
	Dummy = (double*) calloc ( KK, sizeof(double));
	memcpy(Dummy,MeanR,sizeof(double)*KK);
	if (K==1) { iMeanR[0] = 1.0/MeanR[0]; } else { GaussJordan ( Dummy, iMeanR, K);}

	/* Sampled additive effects */
	if (Polygenic)
	{
		sprintf ( Output, "%s_%s_U.txt", Algoname, AnalysisName );
		fp = fopen ( Output, "w" );
		for ( ii=0; ii<Ns; ii++ )
		{
			for ( trait=0; trait<K; trait++ )
				for (additive=0; additive<P; additive++)
					fprintf ( fp, "%f ", Z[trait*P + additive].sampledeffect[ii]);// * Sd[trait]);
			fprintf ( fp, "\n" );
		}
		fclose (fp);
	}

	/* Sampled fixed effects */
	if (F>0)
	{
		sprintf ( Output, "%s_%s_B.txt", Algoname, AnalysisName );
		fp = fopen ( Output, "w" );
		for ( ii=0; ii<Ns; ii++ )
		{
			for (trait=0; trait<K; trait++)
				for (fixed=0; fixed<Ft[trait]; fixed++ )
				{
					if (fixed==0) { temp=Mean[trait]; } else { temp=0.0;}
					fprintf ( fp, "%f ", temp + X[cumFt[trait]+fixed].sampledeffect[ii]);// * Sd[trait]);
				}
			fprintf ( fp, "\n" );
		}
		fclose (fp);
	}

	/* Sampled Lambda */
	if (SEM)
	{
		for (mm=0; mm<M; mm++)
		{
			switch (L[mm].update) {
				case 1:
					sprintf ( Output, "%s_%s_Lambda%d.txt", Algoname, AnalysisName, mm+1 );
					fp = fopen ( Output, "w" );
					for ( ii=0; ii<Ns; ii++ )
						fprintf ( fp, "%f\n", L[mm].sampledlambda [ ii ]); //* Sd[L[mm].down] / Sd[L[mm].up] );
					fclose (fp);
					break;
				case 2:
					sprintf ( Output, "%s_%s_Lambda%d.txt", Algoname, AnalysisName, mm+1 );
					fp = fopen ( Output, "w" );
					for ( ii=0; ii<Ns; ii++ )
					{
						for ( record=0; record<N; record++)
							if (Y[L[mm].up].observations [record] != Missingvalue)
							{
								fprintf ( fp, "%f ", L[mm].sampledlambda [ ii * N + record ]);// * Sd[L[mm].down] / Sd[L[mm].up]);
							}
							else
							{
								fprintf ( fp, "%f ", Missingvalue);
							}
						fprintf ( fp, "\n");
					}
					fclose (fp);
					sprintf ( Output, "%s_%s_Lambda%d_class.txt", Algoname, AnalysisName, mm+1 );
					fp = fopen ( Output, "w" );
					for ( ii=0; ii<Ns; ii++ )
					{
						for ( record=0; record<N; record++)
							if (Y[L[mm].up].observations [record] != Missingvalue)
								fprintf ( fp, "%d ", L[mm].sampledclass [ ii * N + record ] );
							else
								fprintf ( fp, "%d ", -1);
						fprintf ( fp, "\n");
					}
					fclose (fp);
					break;
				case 3:
					sprintf ( Output, "%s_%s_Pi%d.txt", Algoname, AnalysisName, mm+1 );
					fp = fopen ( Output, "w" );
					for(jj=0; jj<Ns; jj++)
					{
						for(ii=0; ii<L[mm].Nsp; ii++)
							fprintf(fp, "%f ", L[mm].sampledPi[jj*L[mm].Nsp+ii]);
						fprintf(fp, "\n");
					}
					fclose (fp);
					sprintf ( Output, "%s_%s_Basis%d.txt", Algoname, AnalysisName, mm+1 );
					fp = fopen ( Output, "w" );
					for(ii=0; ii<L[mm].Nsp; ii++)
					{
						for(record=0; record<N; record++)
							fprintf(fp, "%f ", L[mm].B[record*L[mm].Nsp+ii]);
						fprintf(fp, "\n");
					}
					fclose (fp);
					sprintf ( Output, "%s_%s_Lambda%d.txt", Algoname, AnalysisName, mm+1 );
					fp = fopen ( Output, "w" );
					for(jj=0; jj<Ns; jj++)
					{
						for(record=0; record<N; record++)
						{
							for(ii=0,temp=0.0; ii<L[mm].Nsp; ii++)
								temp += L[mm].sampledPi[jj*L[mm].Nsp+ii]*L[mm].B[record*L[mm].Nsp+ii];
							fprintf(fp, "%f ", temp);
						}
						fprintf(fp, "\n");
					}
					fclose (fp);

					break;
			}
		}
	}

	/* SampledLikelihood */
	sprintf ( Output, "%s_%s_LogLikelihood.txt", Algoname, AnalysisName);
	fp = fopen ( Output, "w" );
	for ( ii=Nswb+1, MeanLoglike=0.0; ii<Ns; ii++)
	{
		for(record=0; record<N; record++)
		{
			temp = log(SampledLikelihood [ii*N+record] );
			fprintf ( fp, "%f ", temp);
			MeanLoglike += temp;
		}
		fprintf (fp, "\n");
	}
	fclose (fp);
	MeanLoglike /= (double)Nsab;

	/* Calculation of DIC and WAIC */
	Products = (double*) calloc ( K, sizeof(double));
	Eri = (double*) calloc ( K, sizeof(double));
	LoglikeFromPostMean = (-0.5*(double)KN)*log(2.0*Pi);
	if (K==1) { LoglikeFromPostMean -= 0.5 * (double) N * log ( MeanR[0]); } else { LoglikeFromPostMean -= 0.5 * (double) N * log(LUdesym (MeanR, Dummy, K));}
	for (record=0; record<N; record++)
	{
		for (trait=0; trait<K; trait++)
			Eri [trait] = Y[trait].meanerrors[record] / (double)Nsab;
		LoglikeFromPostMean += LikelihoodExp(K, iMeanR, Eri, Products);
	}

	PDIC = 2*(LoglikeFromPostMean - MeanLoglike);
	DIC = -2.0*LoglikeFromPostMean + 2.0*PDIC;

	for (record=0, LPPD=0.0, PWAIC2=0.0; record<N; record++)
	{
		for(ii=Nswb+1, temp=0.0; ii<Ns; ii++)
			temp+=SampledLikelihood [ii*N+record];
		temp /= (double)Nsab;
		LPPD += log(temp);

		for(ii=Nswb+1, temp=0.0; ii<Ns; ii++)
			temp+=log(SampledLikelihood [ii*N+record]);
		temp /= (double)Nsab;
		for(ii=Nswb+1, temp2=0.0; ii<Ns; ii++)
			temp2+=pow(log(SampledLikelihood [ii*N+record])-temp,2.0);
		temp2 /= (double)(Nsab-1);
		PWAIC2 += temp2;
	}
	WAIC = -2.0*LPPD + 2.0*PWAIC2;

	sprintf ( Output, "%s_%s_Criterion.txt", Algoname, AnalysisName);
	fp = fopen ( Output, "w");
	fprintf ( fp, "DIC								: %f\n", DIC);
	fprintf ( fp, "WAIC								: %f\n", WAIC);
	fprintf ( fp, "LoglikeFromPostMean				: %f\n", LoglikeFromPostMean);
	fprintf ( fp, "MeanLoglike						: %f\n", MeanLoglike);
	fprintf ( fp, "LogPointwisePredictiveDensity	: %f\n", LPPD);
	fprintf ( fp, "PDIC								: %f\n", PDIC);
	fprintf ( fp, "PWAIC2							: %f\n", PWAIC2);
	fclose (fp);

	/* Free */
	free (Mean); free(Sd);

	if(Polygenic)
	{
		for (additive=0; additive<PK; additive++)
			free(Z[additive].sampledeffect);
		free(SamplediG);
	}

	if (F>0) 
		for (ii=0; ii<F; ii++) 
			free (X[ii].sampledeffect);

	if (SEM)
	{
		for (mm=0; mm<M; mm++)
			switch(L[mm].update)
			{
				case 1:
				free (L[mm].sampledlambda);
					break;
				case 2:
				free (L[mm].sampledlambda);free(L[mm].sampledclass);
					break;
				case 3:
				free (L[mm].sampledPi); free(L[mm].B);
					break;
			}	
	}
	for (trait=0; trait<K; trait++)
		free(Y[trait].meanerrors);

	free(SampledLikelihood); free(Dummy); free(Products); free(Eri);
	free(SamplediR); free(Invertedmatrix); free(Outputmatrix);
	free(MeanR); free(iMeanR);
/*----modify until here for LOO---------*/

	for(trait=0; trait<K; trait++)
		free(Y[trait].observations);
	free (Y);
	if ( Polygenic ) 
	{
		for (additive=0; additive<P; additive++)
			free(Z[additive].covariates); 

		free (Z); free(iA);	free (Va);
	}
	if (CovR==2) free(Ve);
	if (CovR==1) free(ReCor);
	free(Nu); free(S2);
	if (SEM){ free (Structure);	free (L); }
	if (F>0) 
	{ 
		for (ii=0; ii<F; ii++) { free (X[ii].covariates);  }
		free(X); free(cumFt);free(Ft);
	}

	printf ( "Job finished\n");

	return(0);
}
